home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / sspr2.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  7.2 KB  |  233 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE SSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
  5. *     .. Scalar Arguments ..
  6.       REAL               ALPHA
  7.       INTEGER            INCX, INCY, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       REAL               AP( * ), X( * ), Y( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  SSPR2  performs the symmetric rank 2 operation
  17. *
  18. *     A := alpha*x*y' + alpha*y*x' + A,
  19. *
  20. *  where alpha is a scalar, x and y are n element vectors and A is an
  21. *  n by n symmetric matrix, supplied in packed form.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  UPLO   - CHARACTER*1.
  27. *           On entry, UPLO specifies whether the upper or lower
  28. *           triangular part of the matrix A is supplied in the packed
  29. *           array AP as follows:
  30. *
  31. *              UPLO = 'U' or 'u'   The upper triangular part of A is
  32. *                                  supplied in AP.
  33. *
  34. *              UPLO = 'L' or 'l'   The lower triangular part of A is
  35. *                                  supplied in AP.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  N      - INTEGER.
  40. *           On entry, N specifies the order of the matrix A.
  41. *           N must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  ALPHA  - REAL            .
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  X      - REAL             array of dimension at least
  49. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  50. *           Before entry, the incremented array X must contain the n
  51. *           element vector x.
  52. *           Unchanged on exit.
  53. *
  54. *  INCX   - INTEGER.
  55. *           On entry, INCX specifies the increment for the elements of
  56. *           X. INCX must not be zero.
  57. *           Unchanged on exit.
  58. *
  59. *  Y      - REAL             array of dimension at least
  60. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  61. *           Before entry, the incremented array Y must contain the n
  62. *           element vector y.
  63. *           Unchanged on exit.
  64. *
  65. *  INCY   - INTEGER.
  66. *           On entry, INCY specifies the increment for the elements of
  67. *           Y. INCY must not be zero.
  68. *           Unchanged on exit.
  69. *
  70. *  AP     - REAL             array of DIMENSION at least
  71. *           ( ( n*( n + 1 ) )/2 ).
  72. *           Before entry with  UPLO = 'U' or 'u', the array AP must
  73. *           contain the upper triangular part of the symmetric matrix
  74. *           packed sequentially, column by column, so that AP( 1 )
  75. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  76. *           and a( 2, 2 ) respectively, and so on. On exit, the array
  77. *           AP is overwritten by the upper triangular part of the
  78. *           updated matrix.
  79. *           Before entry with UPLO = 'L' or 'l', the array AP must
  80. *           contain the lower triangular part of the symmetric matrix
  81. *           packed sequentially, column by column, so that AP( 1 )
  82. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  83. *           and a( 3, 1 ) respectively, and so on. On exit, the array
  84. *           AP is overwritten by the lower triangular part of the
  85. *           updated matrix.
  86. *
  87. *
  88. *  Level 2 Blas routine.
  89. *
  90. *  -- Written on 22-October-1986.
  91. *     Jack Dongarra, Argonne National Lab.
  92. *     Jeremy Du Croz, Nag Central Office.
  93. *     Sven Hammarling, Nag Central Office.
  94. *     Richard Hanson, Sandia National Labs.
  95. *
  96. *
  97. *     .. Parameters ..
  98.       REAL               ZERO
  99.       PARAMETER        ( ZERO = 0.0E+0 )
  100. *     .. Local Scalars ..
  101.       REAL               TEMP1, TEMP2
  102.       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
  103. *     .. External Functions ..
  104.       LOGICAL            LSAME
  105.       EXTERNAL           LSAME
  106. *     .. External Subroutines ..
  107.       EXTERNAL           XERBLA
  108. *     ..
  109. *     .. Executable Statements ..
  110. *
  111. *     Test the input parameters.
  112. *
  113.       INFO = 0
  114.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  115.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  116.          INFO = 1
  117.       ELSE IF( N.LT.0 )THEN
  118.          INFO = 2
  119.       ELSE IF( INCX.EQ.0 )THEN
  120.          INFO = 5
  121.       ELSE IF( INCY.EQ.0 )THEN
  122.          INFO = 7
  123.       END IF
  124.       IF( INFO.NE.0 )THEN
  125.          CALL XERBLA( 'SSPR2 ', INFO )
  126.          RETURN
  127.       END IF
  128. *
  129. *     Quick return if possible.
  130. *
  131.       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
  132.      $   RETURN
  133. *
  134. *     Set up the start points in X and Y if the increments are not both
  135. *     unity.
  136. *
  137.       IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
  138.          IF( INCX.GT.0 )THEN
  139.             KX = 1
  140.          ELSE
  141.             KX = 1 - ( N - 1 )*INCX
  142.          END IF
  143.          IF( INCY.GT.0 )THEN
  144.             KY = 1
  145.          ELSE
  146.             KY = 1 - ( N - 1 )*INCY
  147.          END IF
  148.          JX = KX
  149.          JY = KY
  150.       END IF
  151. *
  152. *     Start the operations. In this version the elements of the array AP
  153. *     are accessed sequentially with one pass through AP.
  154. *
  155.       KK = 1
  156.       IF( LSAME( UPLO, 'U' ) )THEN
  157. *
  158. *        Form  A  when upper triangle is stored in AP.
  159. *
  160.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  161.             DO 20, J = 1, N
  162.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  163.                   TEMP1 = ALPHA*Y( J )
  164.                   TEMP2 = ALPHA*X( J )
  165.                   K     = KK
  166.                   DO 10, I = 1, J
  167.                      AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
  168.                      K       = K       + 1
  169.    10             CONTINUE
  170.                END IF
  171.                KK = KK + J
  172.    20       CONTINUE
  173.          ELSE
  174.             DO 40, J = 1, N
  175.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  176.                   TEMP1 = ALPHA*Y( JY )
  177.                   TEMP2 = ALPHA*X( JX )
  178.                   IX    = KX
  179.                   IY    = KY
  180.                   DO 30, K = KK, KK + J - 1
  181.                      AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2
  182.                      IX      = IX      + INCX
  183.                      IY      = IY      + INCY
  184.    30             CONTINUE
  185.                END IF
  186.                JX = JX + INCX
  187.                JY = JY + INCY
  188.                KK = KK + J
  189.    40       CONTINUE
  190.          END IF
  191.       ELSE
  192. *
  193. *        Form  A  when lower triangle is stored in AP.
  194. *
  195.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  196.             DO 60, J = 1, N
  197.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  198.                   TEMP1 = ALPHA*Y( J )
  199.                   TEMP2 = ALPHA*X( J )
  200.                   K     = KK
  201.                   DO 50, I = J, N
  202.                      AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
  203.                      K       = K       + 1
  204.    50             CONTINUE
  205.                END IF
  206.                KK = KK + N - J + 1
  207.    60       CONTINUE
  208.          ELSE
  209.             DO 80, J = 1, N
  210.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  211.                   TEMP1 = ALPHA*Y( JY )
  212.                   TEMP2 = ALPHA*X( JX )
  213.                   IX    = JX
  214.                   IY    = JY
  215.                   DO 70, K = KK, KK + N - J
  216.                      AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2
  217.                      IX      = IX      + INCX
  218.                      IY      = IY      + INCY
  219.    70             CONTINUE
  220.                END IF
  221.                JX = JX + INCX
  222.                JY = JY + INCY
  223.                KK = KK + N - J + 1
  224.    80       CONTINUE
  225.          END IF
  226.       END IF
  227. *
  228.       RETURN
  229. *
  230. *     End of SSPR2 .
  231. *
  232.       END
  233.